function [Rall, RA1, RA2, RA12 ] = estim_all_adv_eme(lz, dflz,lx,ldy,lcsx,lcsdy, lcsz, topool,debt_for_ind, debt_for_ind_adv, debt_for_ind_eme,Y, Yadv, Yeme, dX1,dXadv1,dXeme1,indf,indit_f,ind2f,tau_range,S, seladv, seleme )
%INPUTS:
%lx=0;%0; % -1 for no contemporaneous values and no lags of debt and inflation
%ldy=1; % 0 for no lagged dependent variable
%lcsx=-1; %1; %-1; % no contemporaneous values and no lags of cs ave of debt and inflation
%lcsdy=-1; %1; %-1; % no contemporaneous values and no lags of cs ave of dependent variable
%lcsz=-1; %1; %-1; % % no contemporaneous values and no lags of cs ave of indicator variable
%
%topool=false(3,1); % selects which variables to pool in addition to the indicator variable, which is always pooled. See create_regressors_for_threshold_estim for the ordering of individual variables
%debt_for_ind=X(:,:,1);
% tbw


tau_range_adv=tau_range;
tau_range_eme=tau_range;
tic
%% 1
[Regs]=create_regressors_for_threshold_estim(lz, dflz,Y,dX1,lx,ldy, lcsx, lcsdy, indf, lcsz,tau_range, topool, debt_for_ind ); % X(:,:,1) dX1
[phis_AZ1, taus_AZ1]=threshold_estim(Regs, S);
Regsall=Regs;                                         

[Regs]=create_regressors_for_threshold_estim(lz, dflz,Yadv,dXadv1,lx,ldy, lcsx, lcsdy, indf, lcsz,tau_range_adv, topool, debt_for_ind_adv ); % Xadv(:,:,1)
% now replace cs averages with cs averages from the full sample
ip=0;
for i=1:Regsall.N
    if seladv(i,1)
        ip=ip+1;
        if Regsall.posstars{i}(1,1)>0 % cs ave are present
       Regs.x{ip}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2))= Regsall.x{i}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2));
        end
    end
end
[phis_AZ1_adv, taus_AZ1_adv]=threshold_estim(Regs, S);


[Regs]=create_regressors_for_threshold_estim(lz, dflz,Yeme,dXeme1,lx,ldy, lcsx, lcsdy, indf, lcsz,tau_range_eme, topool, debt_for_ind_eme ); % Xeme(:,:,1)
% now replace cs averages with cs averages from the full sample
ip=0;
for i=1:Regsall.N
    if seleme(i,1)
        ip=ip+1;
        if Regsall.posstars{i}(1,1)>0 % cs ave are present
       Regs.x{ip}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2))= Regsall.x{i}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2));
        end
    end
end
[phis_AZ1_eme, taus_AZ1_eme]=threshold_estim(Regs, S);

lra=[ phis_AZ1.lr.xlrmg,phis_AZ1.lr.xlrstd,phis_AZ1.lr.zlrmg,phis_AZ1.lr.zlrstd,phis_AZ1.lr.zxlrmg,phis_AZ1.lr.zxlrstd, phis_AZ1.lr.lambdamg, phis_AZ1.lr.lambdastd];
lradv=[phis_AZ1_adv.lr.xlrmg,phis_AZ1_adv.lr.xlrstd,phis_AZ1_adv.lr.zlrmg,phis_AZ1_adv.lr.zlrstd, phis_AZ1_adv.lr.zxlrmg,phis_AZ1_adv.lr.zxlrstd, phis_AZ1_adv.lr.lambdamg, phis_AZ1_adv.lr.lambdastd];
lreme=[phis_AZ1_eme.lr.xlrmg,phis_AZ1_eme.lr.xlrstd,phis_AZ1_eme.lr.zlrmg,phis_AZ1_eme.lr.zlrstd, phis_AZ1_eme.lr.zxlrmg,phis_AZ1_eme.lr.zxlrstd, phis_AZ1_eme.lr.lambdamg, phis_AZ1_eme.lr.lambdastd];
lr=[lra,lradv,lreme];
klr=size(lr,2);
zl=zeros(1,klr);
RA1=[phis_AZ1.phiall, phis_AZ1.tests.tall, tau_range, phis_AZ1_adv.phiall, phis_AZ1_adv.tests.tall, tau_range, phis_AZ1_eme.phiall, phis_AZ1_eme.tests.tall, tau_range, phis_AZ1.cd,phis_AZ1_adv.cd,phis_AZ1_eme.cd,lr;
    0,0,0, 0,0,0, 0,0,0,0,0,0,0,0,0,zl;
    phis_AZ1.tests.cvals.cvsup,phis_AZ1_adv.tests.cvals.cvsup, phis_AZ1_eme.tests.cvals.cvsup,0,0,0,0,0,0,zl;
    phis_AZ1.tests.cvals.cvave,phis_AZ1_adv.tests.cvals.cvave, phis_AZ1_eme.tests.cvals.cvave,0,0,0,0,0,0,zl; 
    taus_AZ1.tauhat, 0,0, taus_AZ1_adv.tauhat, 0, 0, taus_AZ1_eme.tauhat,0,0,0,0,0,0,0,0,zl ];

%% 2

[Regs]=create_regressors_for_threshold_estim(lz, dflz,Y,dX1,lx,ldy, lcsx, lcsdy, indit_f, lcsz,tau_range, topool, debt_for_ind );
[phis_AZ2, taus_AZ2]=threshold_estim(Regs, S);
Regsall=Regs; 

%debt_for_ind=Xadv(:,:,1);
[Regs]=create_regressors_for_threshold_estim(lz, dflz,Yadv,dXadv1,lx,ldy, lcsx, lcsdy, indit_f, lcsz,tau_range_adv, topool, debt_for_ind_adv );
% now replace cs averages with cs averages from the full sample
ip=0;
for i=1:Regsall.N
    if seladv(i,1)
        ip=ip+1;
        if Regsall.posstars{i}(1,1)>0 % cs ave are present
       Regs.x{ip}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2))= Regsall.x{i}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2));
        end
    end
end
[phis_AZ2_adv, taus_AZ2_adv]=threshold_estim(Regs, S);

%debt_for_ind=Xeme(:,:,1);
[Regs]=create_regressors_for_threshold_estim(lz, dflz,Yeme,dXeme1,lx,ldy, lcsx, lcsdy, indit_f, lcsz,tau_range_eme, topool, debt_for_ind_eme );
% now replace cs averages with cs averages from the full sample
ip=0;
for i=1:Regsall.N
    if seleme(i,1)
        ip=ip+1;
        if Regsall.posstars{i}(1,1)>0 % cs ave are present
       Regs.x{ip}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2))= Regsall.x{i}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2));
        end
    end
end
[phis_AZ2_eme, taus_AZ2_eme]=threshold_estim(Regs, S);

lra2=[phis_AZ2.lr.xlrmg,phis_AZ2.lr.xlrstd,phis_AZ2.lr.zlrmg,phis_AZ2.lr.zlrstd, phis_AZ2.lr.zxlrmg,phis_AZ2.lr.zxlrstd, phis_AZ2.lr.lambdamg, phis_AZ2.lr.lambdastd];
lradv2=[phis_AZ2_adv.lr.xlrmg,phis_AZ2_adv.lr.xlrstd,phis_AZ2_adv.lr.zlrmg,phis_AZ2_adv.lr.zlrstd, phis_AZ2_adv.lr.zxlrmg,phis_AZ2_adv.lr.zxlrstd, phis_AZ2_adv.lr.lambdamg, phis_AZ2_adv.lr.lambdastd];
lreme2=[phis_AZ2_eme.lr.xlrmg,phis_AZ2_eme.lr.xlrstd,phis_AZ2_eme.lr.zlrmg,phis_AZ2_eme.lr.zlrstd, phis_AZ2_eme.lr.zxlrmg,phis_AZ2_eme.lr.zxlrstd, phis_AZ2_eme.lr.lambdamg, phis_AZ2_eme.lr.lambdastd];
lr2=[lra2,lradv2,lreme2];
klr2=size(lr2,2);
zl2=zeros(1,klr2);

RA2=[phis_AZ2.phiall, phis_AZ2.tests.tall, tau_range, phis_AZ2_adv.phiall, phis_AZ2_adv.tests.tall, tau_range, phis_AZ2_eme.phiall, phis_AZ2_eme.tests.tall, tau_range, phis_AZ2.cd,phis_AZ2_adv.cd,phis_AZ2_eme.cd, lr2;
    0,0,0, 0,0,0, 0,0,0,0,0,0,0,0,0, zl2;
    phis_AZ2.tests.cvals.cvsup,phis_AZ2_adv.tests.cvals.cvsup, phis_AZ2_eme.tests.cvals.cvsup,0,0,0,0,0,0, zl2;
    phis_AZ2.tests.cvals.cvave,phis_AZ2_adv.tests.cvals.cvave, phis_AZ2_eme.tests.cvals.cvave,0,0,0,0,0,0, zl2;
    taus_AZ2.tauhat, 0,0, taus_AZ2_adv.tauhat, 0, 0, taus_AZ2_eme.tauhat,0,0,0,0,0,0,0,0, zl2 ];

%% 12

[Regs]=create_regressors_for_threshold_estim(lz, dflz,Y,dX1,lx,ldy, lcsx, lcsdy, ind2f, lcsz,tau_range, topool, debt_for_ind );
[phis_AZ12, taus_AZ12]=threshold_estim(Regs, S);
Regsall=Regs; 

%debt_for_ind=Xadv(:,:,1);
[Regs]=create_regressors_for_threshold_estim(lz, dflz,Yadv,dXadv1,lx,ldy, lcsx, lcsdy, ind2f, lcsz,tau_range_adv, topool, debt_for_ind_adv );
% now replace cs averages with cs averages from the full sample
ip=0;
for i=1:Regsall.N
    if seladv(i,1)
        ip=ip+1;
        if Regsall.posstars{i}(1,1)>0 % cs ave are present
       Regs.x{ip}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2))= Regsall.x{i}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2));
        end
    end
end
[phis_AZ12_adv, taus_AZ12_adv]=threshold_estim(Regs, S);

%debt_for_ind=Xeme(:,:,1);
[Regs]=create_regressors_for_threshold_estim(lz, dflz,Yeme,dXeme1,lx,ldy, lcsx, lcsdy, ind2f, lcsz,tau_range_eme, topool, debt_for_ind_eme );
% now replace cs averages with cs averages from the full sample
ip=0;
for i=1:Regsall.N
    if seleme(i,1)
        ip=ip+1;
        if Regsall.posstars{i}(1,1)>0 % cs ave are present
       Regs.x{ip}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2))= Regsall.x{i}(:,Regsall.posstars{i}(1,1):Regsall.posstars{i}(1,2));
        end
    end
end
[phis_AZ12_eme, taus_AZ12_eme]=threshold_estim(Regs, S);

lra12=[phis_AZ12.lr.xlrmg,phis_AZ12.lr.xlrstd,phis_AZ12.lr.zlrmg,phis_AZ12.lr.zlrstd, phis_AZ12.lr.zxlrmg,phis_AZ12.lr.zxlrstd, phis_AZ12.lr.lambdamg, phis_AZ12.lr.lambdastd];
lradv12=[phis_AZ12_adv.lr.xlrmg,phis_AZ12_adv.lr.xlrstd,phis_AZ12_adv.lr.zlrmg,phis_AZ12_adv.lr.zlrstd, phis_AZ12_adv.lr.zxlrmg,phis_AZ12_adv.lr.zxlrstd, phis_AZ12_adv.lr.lambdamg, phis_AZ12_adv.lr.lambdastd];
lreme12=[phis_AZ12_eme.lr.xlrmg,phis_AZ12_eme.lr.xlrstd,phis_AZ12_eme.lr.zlrmg,phis_AZ12_eme.lr.zlrstd, phis_AZ12_eme.lr.zxlrmg,phis_AZ12_eme.lr.zxlrstd, phis_AZ12_eme.lr.lambdamg, phis_AZ12_eme.lr.lambdastd];
lr12=[lra12,lradv12,lreme12];
klr12=size(lr12,2);
zl12=zeros(1,klr12);


RA12=[phis_AZ12.phiall, phis_AZ12.tests.Wall,  phis_AZ12_adv.phiall, phis_AZ12_adv.tests.Wall,  phis_AZ12_eme.phiall, phis_AZ12_eme.tests.Wall, phis_AZ12.cd,phis_AZ12_adv.cd,phis_AZ12_eme.cd, lr12;
    0,0,0, 0,0,0, 0,0,0,0,0,0,0,0,0, zl12;
    phis_AZ12.tests.cvals.cvsup,phis_AZ12_adv.tests.cvals.cvsup, phis_AZ12_eme.tests.cvals.cvsup,0,0,0,0,0,0, zl12;
    phis_AZ12.tests.cvals.cvave,phis_AZ12_adv.tests.cvals.cvave, phis_AZ12_eme.tests.cvals.cvave,0,0,0,0,0,0, zl12;
    taus_AZ12.tauhat, 0,0, taus_AZ12_adv.tauhat, 0, 0, taus_AZ12_eme.tauhat,0,0,0,0,0,0,0,0, zl12 ];

[s1, s2]=size(RA1);
[a1,a2]=size(RA12);
an=a2-s2;
pom=zeros(7,a2);
Rall=[RA1,zeros(s1,an);
      pom;
      RA2, zeros(s1,an);
      pom;
      RA12];

toc

end

